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We present a nonlinear and non-Markovian random walk model for stochastic movement and 
the spatial aggregation of living organisms that have the ability to sense population density. We 
take into account social crowding effects for which the dispersal rate is a decreasing function of 
the population density and residence time. We perform stochastic simulations of random walk and 
discover the phenomenon of self-organized anomaly (SOA) which leads to a collapse of stationary 
aggregation pattern. This anomalous regime is self-organized and arises without the need for a heavy 
tailed waiting time distribution from the inception. Conditions have been found under which the 
nonlinear random walk evolves into anomalous state when all particles aggregate inside a tiny domain 
(anomalous aggregation). We obtain power-law stationary density-dependent survival function and 
define the critical condition for SOA as the divergence of mean residence time. The role of the initial 
conditions in different SOA scenarios is discussed. We observe phenomenon of transient anomalous 
bi-modal aggregation. 

PACS numbers: 


I. INTRODUCTION 

Aggregation of motile organisms is an example of eco¬ 
logical spatial self-organization due to direct or indirect 
interactions between individuals [I|. Many organisms 
(cells, birds, mammals) have the ability to sense a pop¬ 
ulation density which leads to density-dependent disper¬ 
sal US- This dependence can be explained by var¬ 
ious mechanisms including (i) competition that forces 
an individual to emigrate (positive density-dependence), 
(ii) avoidance of crowded areas by individuals (positive 
density-dependence), (iii) social crowding effects when 
certain areas are attractive to many conspecifics (nega¬ 
tive density-dependence) Q. Density-dependent disper¬ 
sal can be regarded as a behavioral response in which an 
individual changes its rate of jump due to sensing the 
mean population density Si- To model the density- 
dependent dispersal and aggregation, various nonlinear 
local and non-local advection-diffusion equations have 
been used in the literature [M3. Some living organisms 
like amoeboid microorganism Dictyostelium discoideum 
can interact indirectly. They secrete a diffusible attrac- 
tant (signalling molecules) to which individuals respond 
chemotactically. They move towards regions of high con¬ 
centration of attractant and aggregate into a mound. On 
the population level, the standard model for chemical 
interaction of species is a pair of the advection-diffusion- 
reaction equations for the average population density and 
concentration of signalling molecules. There is a huge 
body of literature on chemotactic aggregation modeling 
(see, for example, jlll - IT^ )- 

The microscopic theory of actively moving organisms 
is based on various random walk models |4|-[l4j. One of 
the main characteristics of the continuous time random 
walk (position jump walk) is the escape rate T from the 
point x. For density-dependent dispersal models involv¬ 
ing direct interaction between individuals this rate de¬ 
pends on the population density Q. Indirect interaction 


can be modeled by the dependence of T on concentra¬ 
tion of signalling molecules produced by individuals EO- 
0. An important feature of such random walk models 
is that they are Markovian. However, many transport 
processes are non-Markovian for which the transport op¬ 
erators are non-local in time. Examples include the Levy 
walk that may accelerate aggregation 0 , slow subd- 
iffusive transport that may lead to the phenomenon of 
anomalous aggregation 0 ■ The challenge is to take into 
account both nonlinear density-dependent dispersal and 
non-Markovian anomalous behavior jl9l - f2ll ]. Although 
some research has been done to address the interplay be¬ 
tween nonlinearity and non-Markovian effects |23| . it is 
still an open problem. 

In this paper we consider a nonlinear and non- 
Markovian random walk and propose an alternative 
mechanism of aggregation which we call the regime of 
self-organized anomaly (SOA). It leads to the collapse of 
a standard stationary aggregation pattern and develop¬ 
ment of non-stationary anomalous aggregation. By us¬ 
ing Monte Carlo simulations we show that under certain 
conditions particles performing a non-Markovian random 
walk with crowding effects aggregate inside a tiny do¬ 
main (anomalous aggregation). Important fact is that 
this anomalous regime is self-organized and arises spon¬ 
taneously without the need for a heavy tailed waiting 
time distribution with an infinite mean time from the 
inception. 


II. NONLINEAR AND NON-MARKOVIAN 
RANDOM WALK 

In this section we formulate our nonlinear and non- 
Markovian continuous time random walk model. Instead 
of the waiting time probability density function (PDF) 
we use the escape rate T(r, p) that depends on the resi¬ 
dence time r and the density of particles p. Due to the 
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dependence of the rate T on the residence time r, our 
model is non-Markovian and should involve memory ef¬ 
fects. Our intention is to take into account nonlinear 
social crowding effects and non-Markovian negative ag¬ 
ing. We assume that the random walker has the ability 
to sense the population density p. We model the escape 
rate T as a decreasing function of the density p(x, t) [3,@| 


T(t, p) 


p( t ) 

1 + Ap(x, t) ’ 


(1) 


both sides of the point x as 


e PS(x±a) 

P±( X ) = e /3S(x—a) + e /3S(x+a) ' ( 5 ) 

In this way we introduce the bias of the random walk in 
the direction of the increase of the external field S(x). 
The positive parameter /3 > 0 is the measure of the 
strength of the bias. For small step size a, one can obtain 
the expressions for p± (x) in the continues case 


where A is a positive parameter. This nonlinear function 
describes the phenomenon of conspecific attraction: the 
rate at which individuals emigrate from the point x is re¬ 
duced due to the presence of many conspecifics. This neg¬ 
ative density-dependence can be explained by the various 
benefits of social aggregation like mating, anti-predator 
aggregation, etc. |3| . The rate parameter p (r) is a de¬ 
creasing function of the residence time (negative aging): 


1 B dS . 

p±(x) =-±-—a + o(a). (6) 

Because of the dependence of T on the residence time 
r, it is convenient to define the structured density of par¬ 
ticles £(x, t, i) at time t such that £(x,t, £)AxAt gives 
the number of particles in the space interval (x, x + Ax) 
whose residence time lies in (r, r + At) |23l - t25l |. We con¬ 
sider initial conditions £(x,t, 0) = po(x)6(t), for which 
all particles have zero residence time at t = 0. The total 
density p(x, t) is defined in the standard way 


where po and To are positive parameters. This partic¬ 
ular choice of the rate parameter p (r) has been moti¬ 
vated by non-Markovian crowding: the longer the liv¬ 
ing organisms stay in a particular site, the smaller be¬ 
comes the escape probability to another site. Since 
the escape rate T(r,p) depends on both residence time 
r and t (indirectly through p) we can not define the 
waiting (residence) time PDF. It can be done only for 
the linear case when A = 0. In this case the waiting 
time PDF if) (r) can be defined in the standard way: 

(t) = T(t, 0) exp [— fj T(r, 0)dr] [24[. The particular 
choice © generates the power law distribution: 


V’(r) 


U 0 

P-qTq 

(tq + t) 1+ ^° 


( 3 ) 


For the exponent po < 1, this waiting time probability 
density function has infinite first moment which corre¬ 
sponds to anomalous subdiffusion HIEl. In this paper 
we choose po as 


po > 1 


( 4 ) 


p(x,t) = f €(x,T,t)dT. (7) 
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The balance equation for the density £(x,t, i) for r > 0 
takes the Markovian form 


£(x, r + At, t + At) = £(x, t, t) (1 - T(t, p)At) + o (At), 


where 1 — T(t, p)At is the survival probability during 
At at point x. Since dr/dt = 1, in the limit At —> 0 we 
obtain the following equation 


(K d£ 

dt <9t 


—T(t, p)£. 


( 8 ) 


In what follows we use this equation to determine the 
conditions for the self-orginized anomalous regime. The 
master equation for p can be written as Q 


dp 

dt 


—i(x, t) + p+(x — a)i(x — a, t) 
+P- (x + a)i(x + a,t), 


( 9 ) 


for which the mean waiting time is finite for the linear 
case (A = 0). We do not introduce the anomalous effects 
from the inception as its is done for a classical theory of 
subdiffusive transport fl9i - [2ll ]. 

Regarding the space dynamics, we consider the ran¬ 
dom walk in the stationary external field S(x) on the 
one-dimensional lattice with the step size a. We should 
note that the extension for the two and three dimensions 
is pretty straightforward. When the walker escapes from 
the point x with the rate T(t, p), it jumps to x + a with 
the probability p+(x), and it jumps to x — a, with the 
probability p-(x). For the standard chemotaxis models 
[22| |. it is assumed that the jumping probabilities are de¬ 
termined by the chemoattractant concentration S(x) on 


where i(x,t) = J 0 * T(t, p)£(x, t, t)dr. In general, the ex¬ 
pression for the total escape rate i(x , t) is not known. In 
the Markovian case, when the rate T is independent of 
t, using Eq. 0 we obtain 

i(x,t) = T(p)p(x,t). (10) 

For po < 1 and only for linear case (A = 0), the total 
escape rate takes the form [l^ 

i(x,t) = (T(l-poK°)- 1 D]-^p(x,t), (11) 

where D\~^° is the Riemann-Liouville fractional deriva¬ 
tive. 
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FIG. 1: Population density profiles in linear regime A = 
0 . External field (chemoattractant concentration) is linear: 
S(x) — ax. Parameters are to = 1 , /3 = 1 , a = 0.34 and 
go = 4 . We consider L = 4, and M = 40 , so that a = 0 . 1 . 
The dashed line in the right most panel (coincides with the 
profile) represents the stationary Boltzmann distribution. We 
have used ensemble of 10 6 walkers uniformly distributed at 
time t = 0 on the interval [ 0 . 5 , 2 . 5 ] with zero residence time 
(r = 0 ). The boundaries at x = 0 and x = 4 are assumed to 
be reflecting. 


III. STOCHASTIC SIMULATIONS AND 
SELF-ORGANIZED ANOMALY 

In this section we present the stochastic simulations 
of a nonlinear and non-Markovian continuous time ran¬ 
dom walk along one-dimensional lattice. Because of the 
density-dependent dispersal we intend to model, we can 
not introduce the power law probability density func¬ 
tion for the random residence time as it is done for the 
classical continuous time random walk (CTRW) [ljj-[2l|. 
Therefore we cannot apply the standard Gillespie algo¬ 
rithm f26t . In this paper, we simulate the random walk in 
a domain [0, L] which is divided into M boxes of length 
a = L/M. To calculate the escape rate (JT]) from the box, 
we approximate the density of walkers p in each box as 
the number of walkers in this box divided by the size the 
box. Each walker has its own residence time r which 
determines its escape rate ©■ First, we consider a sta¬ 
tionary linear profile 

S{x) = ax, (12) 

where a is the chemotactic strength parameter. In Fig. [1] 
we present a convergence of the population density profile 
p(x,t) to the stationary distribution p s t(x ) on the inter¬ 
val [0,4] in the linear case {A = 0) for which the walkers 
only sense the external field (concentration of signalling 
molecules) S(x) and not the density of particles p{x,t). 
One can see that the standard aggregation pattern de¬ 
velops. It is easy to show that in the continuous case 
this distribution is the stationary solution of the stan¬ 
dard Fokker-Planck equation: p st (x) ~ exp [ 2(3S(x)]. 

The striking feature of our random walk is that the in¬ 
terplay of nonlinearity © and non-Markovian aging ef¬ 
fect © leads to non-equilibrium phase transition. In the 
nonlinear case (A > 0), when the chemotactic strength 
parameter a is smaller than a critical value a cr , the popu¬ 
lation density convergences to the stationary aggregation 
profile. When the chemotactic strength parameter a is 


greater than a cr , crowding effects induce a collapse of 
the stationary aggregation pattern. The nonlinear ran¬ 
dom walk evolves into the self-organized anomaly. In this 
regime the stationary density profile does not exits and 
all walkers tend to aggregate at the point x m where the 
sub-critical (a < a cr ) stationary density p a t.{x) takes the 
maximum value p s t(x m ). The initial conditions are cho¬ 
sen such that po(x) < p s t{x m ). We discuss the role of 
other initial conditions in section V. Stochastic simula¬ 
tions of the nonlinear and non-Markovian random walk 
are presented in Fig. [2] The top row shows the conver¬ 
gence of the population density to the stationary aggre¬ 
gation profile for a = 0.34. When a < a cr ~ 0.345, we 
observe a development of a stationary aggregation with 
a continues increase in the maximum value of population 
density at the point x m = 4 as the parameter a in¬ 
creases up to a cr . The value of a cr depends on the prop¬ 
erties of the random walk, system size and the boundary 
conditions. Here we do not study this dependence. The 
observed profile is similar to that of Fig. [TJ The only dif¬ 
ference is that the nonlinear crowding effects make the 
value of p s t( 4) greater. However, the drastic change hap¬ 
pens, when the value of a exceeds a cr . The dashed line 
in Fig. [2] represents p st { 4) = 1.5 given by Eq. (l28l) be¬ 
low. Our nonlinear random walk evolves into the non¬ 
stationary SOA which becomes an attractor for random 
dynamics. In this regime all particles eventually concen¬ 
trate in the vicinity of x m = 4. The bottom row in Fig. 
[2] shows this collapse of the density for a = 1. Note that 
similar phenomenon occurs as a result of chemotaxis, the 
so-called chemotactic collapse when all cells aggregate at 
some point. Our explanation of this collapse is different 
from the classical Patlak-Keller-Segel theory in which the 
gro wth of cell density to infinity happens in finite time 

Pill Hi- 

We should note that the self-organized anomaly is an 
universal effect that can occur for any nonuniform exter¬ 
nal field S(x). To demonstrate that the boundary effects 
are irrelevant and to show that the SOA does not de¬ 
pend on the form of S(x), in particular on the derivative 
of S(x) at the point x m , we consider quadratic external 
field with the minimum at the center of the domain [0,4] 

S(x) = -a(x - 2) 2 /2. (13) 

Here a > 0 is the strength parameter. Fig. [3] illustrates 
the phenomenon of the density collapse that takes place 
at the point x m = 2. This shows that SOA is not a bound¬ 
ary effect. 


IV. NONLINEAR MARKOVIAN MODEL 


We should stress the fact that the self-organized 
anomalous regime occurs only for the non-Markovian 
case. In this section we show that there is no anomalous 
collapse for the Markovian nonlinear dynamics when the 
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FIG. 2: Transition to self-organized anomalous regime for the 
linear density of external field S(x) = ax. Transition from 
stationary density to collapsed density occurs when signalling 
strength a exceeds critical value a cr , or in other words when 
the maximum of the stationary density, pst(xm), exceeds the 
critical value p cr = 1.5 given by Eq. (1281) . Here x m = 4 
and po = 4, A = 2. Other parameters, initial and boundary 
conditions are the same as in Fig. [T] The top row shows the 
formation of stationary distribution for a = 0.34 < a cr ■ The 
bottom row illustrates the density collapse that takes place for 
a = 1 > a cr - The critical strength of signalling is estimated 
numerically as a cr — 0.345. 


escape rate T is independent of the residence time r 

Mo 


T(p) = 


r 0 (1 + Ap)' 


(14) 


We start with the Markovian master equation for the 
total density p. Using Eq. ©, we have 

^ = -T{p)p(x,t)+p+(x-a)T(p{x-a,t))p(x-a,t) 
+P-(x + a)T (p( x + a, t)) p(x + a,t ), (15) 


FIG. 3: Same as in Fig. 2 but for quadratic external field 
S( x) = —a(x — 2) 2 /2. The upper row shows the formation of 
stationary distribution for a = 0.64 < cr cr . The bottom row 
illustrates the density collapse that takes place for a = 10 > 
tT cr . Here we consider parameters A = 4 and po = 4. Other 
parameters, initial and boundary conditions are the same as 
in Fig. [2] The critical value of a is estimated numerically to 
be acr — 0.65. Dashed lines represent critical value of the 
density p cr = 3/4 given by Eq. (128II . Notice that in this case 
the collapse takes place at x = 2. 


consider the aggregation process due to negative diffusion 
coefficient [gj. This effect might occur for some negative 
dependence of the escape rate T(p) on the density p. For 
our model with m the last term in Eq. m can be 
modified in a such way that (fTTH) takes the form 


dp 

m 


2/3 


dx 


dS_ 

dx 


D{p)p 


d_ 

dx 


D m (p) 


dp 

dx 


with 


D m (p) 


a 2 Mo 

2r 0 (1 + Ap) 2 


> 0. 


(18) 


where p±(x) is defined by ©• We use the Taylor series 
in m expanding the right hand side in the small a and 
truncate the series at the second term. The equation for 
p takes the form 


= 2/3— 
dt p dx 




+ |^[£(P)P] (16) 


with the nonlinear diffusion coefficient 


Clearly, the modified nonlinear diffusion coefficient 
D m (p) is positive. 


V. THE UNDERLING MECHANISM FOR THE 
SELF-ORGANIZED ANOMALOUS REGIME 


A. Stationary aggregation profile 


D(p) = 


a 2 np) 

2 


q 2 Mo 

2tq (1 + Ap)' 


(17) 


There is no anomalous collapse in this model for any form 
of the signalling concentration S(x). The non-uniform 
stationary solution of (fl6l) with (fI71) gives us a population 
aggregation profile. In the linear case A = 0 we have a 
classical Fokker Planck equation involving the external 
field S(x). We should note that in this paper we do not 


Can we understand the underling mechanism for the 
self-organized anomaly (SOA) observed in our Monte 
Carlo simulations? They show that in the self-organized 
anomalous regime, the standard stationary density pro¬ 
file (aggregation pattern) does not exist. So, it is natural 
first to find a stationary solution £ st (x,T) to the Eq. ® 
from 


d£st 

dr 


—T(r, p st (z))&rf. 













































































































































Using the escape rate © we obtain 


€*t(x,T) = £st(x, 0 ) 


to 


to + T 


Up 

1 + Apst (x) 


This can be rewritten in terms of the density-dependent 
power law survival function 


T,Pst ) = 


TO 




,T 0 +T 

and the stationary arrival rate j s t(x) = £ s t{x, 0) as 


(19) 


Zst(x,T) = j st {x)^(T,p st ). (20) 

In the stationary case the arrival rate of particles j s t{x) 
to the point x and the escape rate of particles from the 
point x are the same: j s t{x) = i s t{x). The stationary 
density p s t{x) can be obtained from ©, (l20l) in the limit 
t —> oo: 


nOO n OO 

Pst{x) = / ^st(x,T)dT = i s t(x) / 'h (t, p st ) dr. 
Jo Jo 


Note that 


( 21 ) 


T ( Pst(x )) 



^(t, p s t)dr 


( 22 ) 


oc=0.5 



-*4- 5 

10 10 10 

t 


FIG. 4: Dependence of the population density at the bound¬ 
ary Xm = 4, p(xm,t), for linear field S(x) = ax. Parameters 
are the same as in Fig. [3 Dashed line mark critical values of 
p cr = 1.5 given by Eq. (1281) . Dashed-dotted curve is the best 
fit for intermediate value of a = 0.4 by logarithmic function 
p(xm,t ) ~ 0.28 ln(t) — 0.22 . 


B. Conditions for self-organized anomaly 


can be interpreted as the expected value of the ran¬ 
dom residence time T whose survival function is given 
by (fT9l) . It should be emphasized that we cannot intro¬ 
duce the power-law survival function as the function of 
non-stationary population density p{x , t). It follows from 
m and ©H) that the stationary escape rate i s t{x) can 
be written in the standard Markovian form 

ist (x) = - 1 .. Pst{x). (23) 

T ( Pst{x)) 


The question arises why the increase of the parameter 
a for the linear external field m or a for the quadratic 
field Eq. (fT3l) above the corresponding critical value a cr 
or a cr leads to a collapse of stationary aggregation pat¬ 
tern (see bottom rows in Figs. [2] and [3]). Our main idea 
is that when a or a gets bigger it leads to the increase 
of the population density at the point x m where density 
p s t{x) takes the maximum value and the divergence of 
the integral 


Using Eq. ©. we write the stationary master equation: 


- i st {x) +p+{x - a)i st (x - a) +p-(x + a)i st (x + a) = 0. 

(24) 

Using probabilities in the limit a —> 0, we obtain 
from (TZ51) and (1M1) the nonlinear stationary advection- 
diffusion equation for the population density p s t{x) : 


2/9 


dx 


S'(x) pst (x) 

d 2 

pst(x) 

_ T(p st (x )) . 

dx 2 

T{p st {x))_ 


= 0. (25) 


It follows from (l25l) that the steady profile p s t(x) 
on the interval [0, L] with the reflecting bound¬ 
aries can be found from the non-linear equation 
Pstix) = N~ lr T (pst) exp [2/3S(x)\, where N is de¬ 
termined by the normalization condition N = 
/ 0 L T (p st (x)) exp [2/3S(x)] dx. This stationary profile 
p s t{x) is illustrated in Fig. [2] for the linear external field 
Eq. (□© and in Fig. [3] for quadratic field Eq. Cl (see 
right most profiles in the top rows). 


'S>(T,p at (x m ))dT. (26) 

In another words: the self-organized anomaly occurs 
when the effective mean residence time T ( x m ) = 

f 0 ^(t, p s t(Xm))dT becomes infinite and the stationary 
Eq. (l?5l) breaks down. The reason why we call this 
regime anomalous is that the divergence of the mean 
waiting time explains anomalous subdiffusive behavior 
of the random walkers [l9l - {2lj . The essential difference 
to the standard CTRW theory is that we use the sta¬ 
tionary density-dependent power law survival function. 
Although SOA is similar to the phenomenon of anoma¬ 
lous aggregation []j§] or the accumulation of subdiffusive 
particles in one of two infinite domains with two different 
values of anomalous exponents [27|, it is essentially dif¬ 
ferent. Anomalous conditions are not imposed by power 
law waiting time PDF © with the anomalous exponent 
po < I- Anomalous regime is self-organized for po > 1 
through the nonlinear interactions of random walkers due 
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FIG. 5: The role of the initial conditions for the quadratic 
external field S(x) with different field strength a. We have 
considered parameters to = 1, (5 = 1, A = 4, po = 4, L = 
4 and M = 400, so that a = 0.01. In contrast to Fig. 0 
the initial distribution of particles is chosen to violate the 
condition in Eq. m • We have used ensemble of 10 6 walkers 
uniformly distributed at time t = 0 on the interval [0.5, 0.6] 
with zero residence time (r = 0). The boundaries at x = 0 
and x = 4 are assumed to be reflecting. The dashed lines 
correspond to the critical value of the density p cr = 3/4, see 
Eq. (1281) . Curve 1 corresponds to t = 10 4 and curve 2 to 
t = 7 x 10 4 . 


to social crowding effects described by ©• 

Substitution of the survival function © into © 
gives 


t 0 [1 + Ap st {x m )] 
Mo - 1 - Ap s t(x m )' 


(27) 


The divergence of T(x m ) gives the critical density p cr : 


Per — 



(28) 


One can also write the critical condition as 7 = 1, where 


7 = 


Mo 


1 4“ Apst (Xm ) 


Mo > 1- 


(29) 


In particular, for the linear external field S(x) = ax in 
the interval [ 0 ,4] the stationary density p s t{x) has a max¬ 
imum value at the point x m = 4. We can find the critical 
value a cr as 


lim 

Oi—tOicr 


Po 

1 T Ap st (x m ) 


= 1 , 


Mo > 1- 


(30) 


Numerical simulations presented in Fig. [2] and [3] are in 
excellent agreement with Eq. ( 1251 ) . In Fig. 0 ] we illustrate 
the time evolution of the population density p(x m ,t) at 
the boundary x m = 4 for linear field S(x) with different 
strength of the signalling a. Transition to self-organized 
anomalous regime is observed by the transition of the 
density through the critical value given by Eq. (l28l) . For 
intermediate values of a, we find that the density grows 
as In t see Fig. 0] 


C. The role of the initial conditions. 

In our numerical simulations we have used the initial 
conditions for which all walkers have zero residence time 


(no aging effects) and the maximum value of the initial 
density po(x) obeys the inequality 

max p 0 (x) < p cr , (31) 

where p cr is given by Eq. (1281) . However, when this con¬ 
dition is violated the SO A regime could have different 
scenarios. In Fig. [5] some preliminary results of numer¬ 
ical simulations are shown. We observe three scenarios 
depending on the strength of the external field er. The 
left most panel shows the anomalous aggregation of the 
walkers in the region of non-zero initial density [0.5, 0.6]. 
The peak of the density in this region is infinitely increas¬ 
ing (this increase is not captured in Fig. [5], only the peak 
is visible). The density in the other region reaches a quasi 
stationary form approximately represented by curve 2 for 
t = 7 x 10 4 . Note, however, that it is not a true station¬ 
ary state since the particles continue to accumulate in the 
interval [0.5,0.6]. The middle panel demonstrates simi¬ 
lar scenario for a = 1. However, now the field is strong 
enough to temporarily accumulate particles at the crit¬ 
ical pointa; m = 2. Our assumption is that similarly to 
the previous case, in the long time limit these particle 
will be concentrated in the region of non-zero initial con¬ 
ditions. A different scenario occur for large strength of 
the external field. The right most panel shows this sit¬ 
uation for a = 5. In this case the external field is able 
to strongly attract particles to its critical point x m = 2 . 
At some time the density at this point exceed the criti¬ 
cal value Eq. ( 1251 ) . Therefore, in this case the anomalous 
aggregation of the density is observed at two places, at 
the region of non-zero initial conditions and at the criti¬ 
cal point x m = 2. We call this phenomenon as transient 
anomalous bi-modal aggregation. 


VI. CONCLUSIONS 

We have discovered the phenomenon of self-organized 
anomaly (SOA) which takes place in a population of par¬ 
ticles performing nonlinear and non-Markovian random 
walk. This random walk involves social crowding effects 
for which the dispersal rate of particles is a decreasing 
function of the population density and residence time 
II- Monte Carlo simulations shows that the regime of 
self-organized anomaly leads to a collapse of a station¬ 
ary aggregation pattern when all particles concentrate 
inside a tiny region of space and form a non-stationary 
high density cluster. The maximum population density 
slowly increases with time as lnt. We should note that 
the anomalous regime is self-organized and arises spon¬ 
taneously without the need for introduction of the power 
law waiting time distribution with infinite mean time. 
Only in a stationary case one can obtain a power-law 
density-dependent survival function and define the crit¬ 
ical condition as the divergence of mean residence time. 
SOA gives a new possible mechanism for chemotactic col¬ 
lapse in a population of living organisms as an alternative 
to the celebrated Patlak-Keller-Segel theory PEIS. 
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The crossover from the standard stationary aggregation 
pattern to a non-stationary anomalous aggregation as the 
strength of chemotactic force increases can be interpreted 
as non-equilibrium phase transition. Our theory can 
be used to explain various anomalous aggregation phe¬ 
nomenon including accumulation of phagotrophic pro- 
tists in attractive patches where they become almost im¬ 
mobile [28]. It would be interesting (i) to analyze ex¬ 


tensions of our model by taking into account additional 
effects like volume feeling preventing unlimited density 
growth t 23i and (ii) apply our model for the analysis of 
how self-organization and an anomalous cooperative ef¬ 
fect arise in social systems [29}. 

This work was funded by EPSRC grant EP/J019526/1. 
The authors thank Steve Falconer for the helpful sugges¬ 
tions. 
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